###   Set the root directory of the project here, eh?
setwd("/Users/michael/Dropbox/carbonefficiency")

###########################
##   specify these please

m <- 2         ## number of inputs
s <- 2         ## number of good outputs
h <- 1         ## number of bad outputs
nd <- 2        ## number of nondiscretionary variables  

##  Thank you... dope


echo=FALSE;

#library(foreign);
library(boot);
library(MASS)
#library(np)
#library(nonparaeff)
#library(Benchmarking)
library(lpSolve)
library(Hmisc)
library(car)


###   work station
dat.mat <- read.table((paste(getwd(),"/data/final/final_data2011.csv",sep="")), header = TRUE, sep = ",", stringsAsFactors = FALSE)

##  run this part
############################
##  subset only variables needed
##  data table needs to be in order of inputs, good outputs, bad outputs, non-discretionary vars
##  with each column being one input or output


attach(dat.mat)
dat.tab <- data.frame(instnm, CarnClass,    ##  identifiers
                      InstrExp = InstrExp/1000, ResExp = ResExp/1000,     ##  inputs
                      FTE, TotRes = TotRes/1000,      ##  good outputs
                      totemiss12,           ##  bad outputs
                      hdd, cdd   ##  non-discretionary vars
                      )       
detach(dat.mat)
dat.mat <- dat.tab


res.schools <- subset(dat.mat, (CarnClass == "R1" | CarnClass == "R2" | CarnClass == "R"))
R1.schools <- subset(dat.mat, (CarnClass == "R1" ))
R2.schools <- subset(dat.mat, (CarnClass == "R2" ))
AS.schools <- subset(dat.mat, (CarnClass == "AS" ))


var.sum.table <- function(variable){
	all <- get(paste("dat.tab$",variable, sep = ""))
	data.mat <- data.mat[,3:ncol(data.mat)]
	sum.df <- data.frame(variable = colnames(data.mat),
	                mean = colMeans(data.mat),
	                std.dev = apply(data.mat,2,sd),
	                min = apply(data.mat,2,min),
	                max = apply(data.mat,2,max) 
	                 )
	                 print(sum.df)
	                 apply(data.mat, 2, hist)
	sumtex <- latex(format.df(sum.df,dec = 2), file = "", rowname = NULL)
	return(sumtex)
}

sum.table <- function(data.mat){
	data.mat <- data.mat[,3:ncol(data.mat)]
	sum.df <- data.frame(variable = colnames(data.mat),
	                mean = colMeans(data.mat),
	                std.dev = apply(data.mat,2,sd),
	                min = apply(data.mat,2,min),
	                max = apply(data.mat,2,max) 
	                 )
	                 print(sum.df)
	                 apply(data.mat, 2, hist)
	sumtex <- latex(format.df(sum.df,dec = 2), file = "", rowname = NULL)
	return(sumtex)
}



####   Nothing to see here... move along
####   These aren't the droids you're looking for
###################################
###   plot Instruction expenditure

# create a vector to store the variable names
var.names <- c("All schools", "LibArts", "R1", "R2", "Research") #,"HC1","HC2","HC3","HC4","HC4m","HC5")

# add the data
cent <- rep(mean(dat.tab$InstrExp),5)        # center for all plotted lines
est[2] <- coef(model.fgls)[11]
est[3] <- coef(model.npgls)[11]
#est <- c(est.npgls, est.fgls, est)
se <- c(se.ols[11], se.fgls[11], se.npgls[11], se.hc0[11],se.hc1[11],se.hc2[11],
		se.hc3[11],se.hc4[11],se.hc4m[11],
		se.hc5[11])           			# Store the std. errors 


pdf(file="CoefficientPlot.pdf")



# set the graphical parameters
par(
  family = "serif",  # I don't plot in anything but serif
  oma = c(0,0,0,0),  # Since it is a single plot, I set the outer margins to zero.
  mar = c(5,10,4,2)  # Inner margins are set through a little trial and error.
  )

# create an empty plot for total customization
plot(NULL,                              # create empty plot
  xlim = c(-0.5, 0.10),                  # set xlim by guessing
  ylim = c(.7, length(var.names) + .3), # set ylim by the number of variables
  axes = F, xlab = NA, ylab = NA)       # turn off axes and labels

# loop over a counter the length of the estimate vector

for (i in 1:length(est)) {

	points(est[i], i, pch = 19, cex = .95)                   # add the points to the plot
	
	lines(c(est[i] + qnorm(0.975)*se[i], 
			est[i] - qnorm(0.975)*se[i]), c(i, i))         	# add the 95% confidence intervals
	
	lines(c(est[i] + qnorm(0.95)*se[i], 
			est[i] - qnorm(0.95)*se[i]), c(i, i), lwd = 3)  # add the 90% confidence intervals
	
	text(-0.6, i, var.names[i], xpd = T, cex = .8)          # add the variable names

}

# add axes and labels
axis(side = 1)                                                         # add bottom axis
abline(v = 0, lty = 3, col = "grey")                                   # add verticle line
abline(v=est[1],lty=4,col="grey")
#mtext(side = 1, "Heteroscedasticity Robust Confidence Intervals", line = 3)  # label bottom axis
mtext(side = 3, "90 and 95% Confidence Intervals Across \n Competing Estimators for Easiness", line = 1)   # add title
box()                                                                                                   # add lines around the plot
dev.off()




